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1 . The  Operational  Approach 

There  are  essentially  two  possible  approaches  to  physical 
geodesy  (as  also  to  other  natural  sciences):  they  might  be  called 
the  model  approach  and  the  operational  approach.  Essentially,  the 
first  approach  starts  from  a theory,  the  second  from  the  observa- 
tions. Obviously,  the  two  approaches  are  closely  related  to  the 
deductive  method  and  the  inductive  method  in  the  natural  sciences. 

In  the  model  approach,  one  starts  from  a mathematical  mo- 
del or  from  a theory  and  then  tries  to  fit  this  model  to  reality, 
for  instance  by  determining  the  parameters  of  this  model  from  ob- 
servations, The  classical  geodetic  example  are  the  centuries-old 
attempts  to  determine  the  parameters  of  an  earth  ellipsoid  by 
observation,  from  the  old  grade  measurements  to  modern  satellite 
observations . 

Perhaps  the  most  elaborate  form  of  this  model  approach  is 
the  boundary-val ue  problem  of  physical  geodesy  in  the  formulation 
of  Molodensky.  It  has  a mathematically  enormously  interesting  and 
deep  theory  and  is  practically  highly  significant,  as  the  many 
gravimetric  geoid  determinations  and  computations  of  deflections 
of  the  vertical  show.  However,  this  approach  has  its  weeknesses: 
the  required  continuous  gravity  coverage  is  practically  not 
realizable;  on  the  other  hand,  many  other  important  data  cannot 
be  incorporated  into  this  theory.  The  model  selects  its  data. 

At  present  we  have  a great  number  of  geodetic  measurements 
of  very  different  types,  from  terrestrial  angle  and  distance 
measurements  to  satellite  data  of  various  kinds.  The  question 
arises:  how  can  we  use  and  combine  all  these  data  in  the  best 
possible  way.  This  is  the  operational  approach. 

Let  us  summarize.  In  the  model  approach  one  asks:  how  can 
I best  determine  my  model  by  suitable  observations?  In  the  opera- 
tional approach  one  asks:  how  can  I make  best  use  of  all  my  ob- 
servations? 

As  a matter  of  fact,  the  two  approaches  do  not  compete 
with  each  other;  each  one  incorporates  important  aspects,  and  the 


2 


two  approaches  mutually  complement  each  other. 

The  operational  approach  to  physical  geodesy  has  come  up 
at  a relatively  recent  date,  when  a huge  number  of  measurements 
of  new  types  was  available  and  when  it  turned  out  that  the  clas- 
sical, especially  the  gravimetric  approach  failed  to  give  a com- 
plete answer  in  view  of  the  lack  in  gravity  data. 

In  geometrical  geodesy  already  least-squares  adjustment 
is  in  the  spirit  of  an  operational  approach  (how  can  I best  use 
all  my  measurements).  In  physical  geodesy,  operational  methods 
have  been  known  under  the  names  "least-squares  collocation", 
"integrated  geodesy",  "operational  geodesy".  All  these  methods 
are  very  similar;  they  all  aim  at  an  adequate  treatment  of  the 
gravity  field,  in  addition  to  an  adjustment  of  measuring  errors. 
They  all  use  quadratic  minimum  principles  incorporating  not  only 
the  measuring  errors,  but  also  the  anomalous  gravity  field. 

The  purpose  of  the  present  report  is  a systematic  treat- 
ment of  the  operational  approach.  Obviously,  the  determination 
of  the  gravity  field,  of  station  coordinates,  etc.,  from  discrete 
measurements  does  not  in  general  admit  a unique  solution;  it  is 
an  "improperly  posed  problem".  The  question  is:  Which  practically 
feasible  possibilities  exist?  Is  there  a genuine  alternative  to 
least-squares  collocation  and  related  methods,  which  is,  perhaps, 
even  simpler? 

This  question  is  treated  by  using  modern  techniques  for 
solving  improperly  posed  problems.  They  all  seem  to  converge  on 
collocation  by  means  of  kernel  functions  and  least-squares  collo- 
cation; but  we  shall  also  try  to  point  out  alternatives. 


2.  Measurements  as  Nonlinear  Functionals 

Every  geodetic  measurement  depends: 

1.  on  one  or  several  points  In  space; 

2.  on  the  earth’s  gravitational  field. 
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Symbolically  we  may  write: 

1 = F (X , V ) . (2-1) 

Here  1 denotes  the  measurement  under  consideration,  V denotes 
the  gravitational  potential,  and  the  vector  X comprises  the 
coordinates  of  the  points  to  which  the  measurement  refers,  and 
possibly  other  parameters.  For  instance,  if  we  have  two  points 
P and  Q and  if  we  use  rectangular  coordinates  xyz  referred 
to  some  Cartesian  reference  system,  then 


T denoting  the  transpose  (in  the  sequel,  vectors  will  be  column 
vectors  unless  the  contrary  is  stated). 

The  symbol  F denotes  any  functional  dependence  on  X 
and  V.  With  respect  to  X,  it  is  an  ordinary  function;  but  it 
is  not  necessarily  an  ordinary  function  of  V but  may  involve 
first  and  higher  derivatives  of  V,  integrals,  etc.  In  the  ter- 
minology of  functional  analysis,  F is  a (nonlinear)  functional 
of  X and  V. 

Denote  the  number  of  components  of  x by  m;  then  X 
may  be  said  to  belong  to  m -dimensional  Euclidean  space  R . The 

m 

function  V may  be  considered  to  belong  to  some  set,  or  space, 

H of  harmonic  functions.  Then,  in  the  jargon  of  modern  mathema- 
tics, the  functional  F is  a mapping  of  the  product  space 
R x H into  R,  the  real  number  line: 

m 

F:  Rm  x H - R , (2-3) 

which,  in  the  language  of  plain  mortals,  means  simply  that  F 

associates,  to  each  harmonic  function  from  the  set  H and  to 

each  vector  R , a real  number  which  represents  the  numerical 
m 

value  of  the  observation  1. 


I 

\ 

i 

i 
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More  intuitively  we  may  say  that  F is  nothing  else  but 
a prescription  for  computing  a number  1 from  a given  vector  X 
and  a given  function  V:  if  X and  V are  supposed  to  be  known, 
then  it  must  be  possible  to  find,  in  an  unambiguous  way,  the  va- 
lue of  1 . In  other  terms,  F denotes  an  operation  to  be  per- 
formed on  X and  V,  the  result  of  which  is  a real  number. 

Thus,  functionals  are  special  cases  of  operators:  for  an 
operator,  the  result  of  the  operation  may  be  quite  general:  a 
vector  or  another  function;  if  the  result  is  a real  number,  the 
operator  is  called  a functional.  Best  known  from  functional  ana- 
lysis (cf.  Kantorovich  and  Akilov,  1964;  Meissl,  1975;  Tscherning, 
1978)  are  linear  operators  and  linear  functionals.  Nonlinear 
operations  are  basic  in  the  so-called  "Modern  Analysis"  (Dieu- 
donne,  1960;  Loomis  and  Sternberg,  1968). 

Our  functionals  (2-1)  will,  in  general,  be  nonlinear;  in 
the  next  section,  we  shall  describe  how  they  can  be  linearized. 

The  physical  meaning  of  such  functionals  F will  be  clear 
from  the  examples  given  below. 

Instead  of  the  gravitational  potential  V,  we  may  also  use 
the  gravity  potential  W,  defined  in  the  usual  way  by 

W * V + \ to2 (x2  + y2)  , (2-4) 

a)  denoting  the  angular  velocity  of  the  earth's  rotation  and  the 
z -axis  coinciding  with  the  earth's  axis  of  rotation.  For  the 
present  purpose  it  is  appropriate  to  regard  w as  constant  and 
the  rotation  axis  as  Invariable  with  respect  to  the  earth's  body: 
the  very  small  deviations  from  this  Idealized  situation  can  easi- 
ly be  taken  into  account  by  corrections,  without  changing  the 
essential  picture. 

Throughout  this  report,  the  system  xyz  will  thus  be  de- 
fined as  follows:  the  origin  Is  at  the  earth's  center  of  mass, 
the  z -axis  coincides  with  the  (mean)  rotation  axis,  and  the  x- 
axis  goes  through  the  (mean)  Greenwich  meridian. 
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Then,  instead  of  (2-1),  we  have 

1 = F( X , W) . (2-5) 

As  W is  expressed  in  terms  of  V and  X by  (2-4),  the  rela- 
tions (2-1)  and  (2-5)  are  equivalent;  as  a matter  of  fact,  the 
letter  F denotes  different  functionals  in  each  of  the  two  ca- 
ses . 

Let  us  now  illustrate  these  abstract  considerations  by 
means  of  concrete  examples. 

Astronomical  and  gravimetric  observations.  The  gravity 
vector  G is  expressed  in  terms  of  gravity  g and  astronomical 
latitude  $ and  longitude  A by  means  of  the  well-known  rela- 
tion 


- € = 


gcos$cosi\1 
gcos*s i nA 
[_gs  i n* 


(2-6) 


(cf.  Heiskanen  and  Moritz,  1967,  p . 5 7 ) . On  the  other  hand,  G 
is  the  gradient  of  the  gravity  potential  W : 


G = grad  W = 


W 


(2-7) 


denoting  the  partial  derivative 


W = | « 

x 3 X 


(2-8) 


and  similarly  for  Wy  and  Wz.  Comparing  (2-6)  and  (2-7)  and 
solving  for  *,  A,  and  g we  obtain 


* * tan 


•l  -W, 


/w2  + W2 

x y 


(2-9) 


isJfcuS)  JfcfeL 


6 


A 


tan 


X 


(2-10) 


g = /w2  + W2  + W2 
3 x y z 


(2-11) 


These  equations  have  the  form  (2-5):  they  express  the  observables 
$ , A,  g in  terms  of  the  potential  W , not  as  ordinary  functions 
of  W , but  as  nonlinear  functionals  involving  the  operation  of 
differentiation.  Let  X denote  the  coordinate  vector  of  the  ob- 
servation station: 


(2-12) 


Then,  as  , Wy  , Wz  are  functions  of  x,  y,  z,  the  expressions 
(2-9)  to  (2-11)  do,  in  fact,  also  depend  on  X , in  agreement  with 
(2-5). 

Angle  and  distance  measurements.  The  observables:  azimuth 
A , zenith  distance  Z , and  distance  S between  two  points  P 
and  Q , can  be  expressed  in  terms  of  the  coordinate  differences 


AX  = XQ  - xp  , 

Ay  = yQ  - yp  , 


as  follows  (Heiskanen  and  Moritz,  1967,  p.219): 


. „ t -l  - AxsinA  + AycosA 

- Axs1n*cosA  - Aysin4>sinA  + azcos$ 

Z * co--i  AxcosftcosA  + AycosftsinA  + Azsins 

/ax2  + Ay2  + az2 

s * /ax2  + Ay2  + az2 


(2-13) 


(2-14) 

(2-15) 

(2-16) 


1 


I 

I 

1 


f 

f 

Again,  these  equations  have  the  form  (2-5);  the  vector  X is 
now 


comprising  the  coordinates  of  both  points  P and  Q,  and  the 
dependence  on  the  potential  W is  implicitly  through  <t>  and  A 
as  expressed  by  (2-9)  and  (2-10);  hence  A and  Z are,  in  fact, 
nonlinear  functionals  of  W.  Note  that  these  observables  depend 
on  the  target  point  Q only  because  its  coordinates  enter  into 
Ax,  Ay,  Az;  on  the  observation  station  P they  depend  in  the 
same  way,  but  there  is  an  additional  dependence  on  P because 
i>  and  A,  and  hence  Wx,  Wy,  Wz,  refer  to  this  point. 

A measured  horizontal  angle  u>  may  be  considered  as  the 
difference  between  two  azimuths: 


u>  = A2  - A,  , (2-18) 

measured  at  an  observation  station  P to  two  targets  and  Q2< 

Both  azimuths  Aj  and  A2  may  be  expressed  by  (2-14);  the  resul- 
ting expression  for  u>  clearly  involves  the  coordinates  of  P, 

Qj,  and  Q2,  so  that,  in  the  present  case,  the  vector  X consists 
of  9 components,  which  are  the  coordinates  of  these  three  points; 
we  again  get  a nonlinear  functional  of  form  (2-5). 

It  goes  without  saying  that  the  functional  expression  t 

(2-16)  for  the  distance  (2-16)  Is  also  a special  case  of  (2-5), 
in  which  there  simply  is  no  factual  dependence  on  W;  measured 
straight  distances  between  two  fixed  points  do  not  depend  on  the 
gravity  field. 

Satellite  observations.  Consider  a distance  S measured 
from  a ground  station  P to  a satellite  Q by  laser  or  radar. 

(Cf.  Fig.  2-1;  for  a non-technl cal  and  compact  review  of  various 
techniques  see  (Cordova,  1977).) 


e direction  observation 

fj 

S distance  measurement 

dS/dt  doppler  observation 
h satellite  altimetry 

ds/dt  sate  1 1 i te-to-satel 1 1 te  tracking  (doppler) 

* 1 


Figure  2-1.  Satellite  techniques 
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Such  a distance  can  again  be  represented  by  (2-16)  but,  if  we 
operate  in  the  orbital  mode,  the  coordinates  of  Q can  be  fur- 
ther expressed  by  the  six  orbital  elements  p , p , ....  p of 
some  reference  orbit  and  the  coefficients  J and  K of  the 

nm  run 

expansion  of  the  earth's  gravitational  potential  V in  terms  of 
spherical  harmonics.  Thus  S will  have  the  form  of  some  function 


S = S(xp,yp,zp;p1,p2 


,pc;J  , K ) 

r 6 nm  nm ' 


(2-19) 


This  is  a functional  of  form  (2-1).  The  vector  X is  given  now 
by 


X 


VYp’VPi’Ps’ 


(2-20) 


it  comprises  station  coordinates  and  orbital  parameters.  The  sphe- 

rical -harmonic  coefficients  J and  K may  be  expressed  in 

nm  nm  r 

terms  of  V by  well-known  integral  formulas  (of  type  of  eq.(l-76) 
of  (Heiskanen  and  Moritz,  1967,  p.31)),  which  explains  the  func- 
tional dependence  on  V. 

The  change  of  distance  S with  respect  to  time  t,  that 
is,  the  range-rate  d3/dt,  can  be  measured  by  doppler  observations. 
By  integrating  dS/dt  with  respect  to  t from  t1  to  t2,  one 
obtains  distance  differences  S2  - S . By  photographing  the  sa- 
tellite against  the  background  of  stars  one  finds  the  right  as- 
cension and  the  declination  of  the  spatial  direction  PQ,  or  in 
other  terms,  the  unit  vector  e of  this  direction  (Fig.  2-1).  All 
these  observables  have  the  same  mathematical  structure  as  (2-19): 
they  are  again  functionals  of  type  (2-1),  the  vector  X being 
given  by  (2-20). 

Satellite  altimetry  can  be  considered  to  measure  the  height 
h of  a satellite  above  the  geoid:  the  ocean  surface  reflects  a 
radar  signal  emitted  by  the  satellite,  and  under  idealized  condi- 
tions, this  surface  coincides  with  the  geoid.  We  claim  that  h 
can  again  be  expressed  as  a functional  of  type  (2-1),  with 
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^ * ^Pl»P2****»P 6^  * (2-21) 

This  is  true  if  it  is  possible,  given  X and  the  potential  func- 
tion V(x,y,z),  to  compute  h.  In  fact,  assume  the  gravitational 
potential  V(x,y,z)  to  be  known  as  a function  of  position  at  all 
points  outside  and  on  the  earth's  surface.  Then  the  gravity  po- 
tential function  W is  also  known  by  (2-4),  and  consequently  the 
geoid  is  an  equi potenti al  surface 

W(x,y,z)  = WQ  = const.  (2-22) 

Now,  the  satellite  orbit  can  be  computed  from  the  parameters  pk 
of  the  reference  orbit  and  the  gravitational  potential  V,  and 
the  position  Q of  the  satellite  along  the  orbit  is  uniquely  de- 
termined by  giving  the  corresponding  instant  t (which  we  assume 
to  be  known).  Thus  both  the  geoid  and  the  satellite  position  Q 
are  determined,  and  so  is  h,  as  the  length  of  the  perpendicular 
from  Q to  the  geoid.  Therefore,  by  the  definition  of  the  func- 
tional (2-1)  given  above,  the  satellite  altimeter  measurement  h 
is,  in  fact,  such  a functional. 

The  data  of  satell 1 te-to-satel li te  tracking  are  time  chan- 
ges of  the  distance  s between  two  satellites  (Fig.  2-1).  Such 
a range  rate  ds/dt  Is  again  measured  by  means  of  the  doppler 
principle.  At  present  one  generally  uses  one  high  and  one  low  sa- 
tellite, but  the  use  of  two  low  satellites  which  are  close  to 
each  other  Is  also  possible.  Considerations  analogous  to  the  pre- 
ceding ones  make  It  obvious  that  ds/dt  has  again  the  form  (2-1), 
the  vector  X comprising  now  the  6+6  elements  of  the  two  re- 
ference orbits. 


Satellite  gradlometry  Is  designed  to  measure  elements  (or 
linear  combinations  of  elements)  of  the  second-order  gradient 
tensor 
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V 

V 

XX 

xy 

xz 

V 

V 

xy 

yy 

yz 

V 

V 

xz 

yz 

zz 

(2-23) 


which  is  a symmetric  matrix  formed  by  the  second  derivatives  of 
the  potential  V with  respect  to  the  coordinates  xyz.  Any  se- 
cond-order gradient,  say  V , depends  on  position: 

xz 


V 

xz 


V ( x ,y , z ) . 

xz 


It  has,  therefore,  the  form  (2-1),  with 


(2-24) 


(2-25) 


the  prescription  for  computing  the  functional  F in  the  present 

case  consists  in  differentiating  V with  respect  to  x and  z 

and  taking  V at  the  point  with  coordinates  (2-25). 
xz 

Very -long- base! i ne  interferometry  measures  the  delay  t, 
with  which  a radio  signal  emitted  from  an  extragal acti c radio 
source  is  received  at  two  different  places  P and  Q (Macdoran, 
1973).  By  multiplying  T with  the  light  velocity  c we  get  the 
projecti on 


D 


PQ.e 


(2-26) 


of  the  vector  PQ  connecting  the  two  points  onto  the  direction 
(supposed  known)  to  the  radio  source  represented  by  the  unit  vec- 
tor e.  Similarly  to  a distance  measurement  (2-16),  D does  not 
depend  on  the  gravity  field,  and  we  have  a special  case  of  (2-1) 
with  X being  given  by  (2-17)  and  with  no  explicit  dependence 
on  V . 

These  examples  should  make  it  obvious  that  al 1 geodetic 
measurements,  without  exception,  can  be  represented  as  functionals 
(2-1)  or  (2-5).  This  simple  and  general  fact  will  be  basic  for  the 
considerations  to  follow. 
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It  is  clear  that  we  have  taken  into  account  only  the  geo- 
metrical and  gravitational  structure  of  the  problem.  We  have  ab- 
stracted from  random  and  systematic  errors,  nongravi tational 
effects,  etc.  Random  errors  will  be  considered  later  in  this  re- 
port, and  systematic  effects  are  assumed  to  have  been  removed  by 
appropriate  corrections.  If  necessary,  systematic  parameters  can 
be  included  in  the  vector  X in  (2-1)  or  (2-5). 


3.  Linearization 


Every  observation  1 gives  an  equation  of  type  (2-1)  or 
(2-5).  We  thus  obtain  a system  of  functional  equations 


\ = Ft(X,W)  , 


12  » F2(X,W)  , 


(3-1) 


lq  = Fq(X,W)  , 


which  are  to  be  solved  for  the  unknown  parameters 
tential  function  W. 


X and  the  po- 


Since  the  functionals 


Fl’  F 2 * 


are  non-1 inear, 


the  system  (3-1)  is  very  difficult  to  handle  directly.  The  usual 
procedure  with  difficult  nonlinear  problems  is  to  linearize  them 
by  Taylor's  theorem. 


Let  us  introduce  an  approximate  value  XQ  for  the  vector 


X and  an  approximation  U to  the  gravity  potential  W.  The  func- 
tion U Is  called  the  normal  potential  ; it  is  generally  taken  to 
be  the  external  gravity  potential  of  an  equipotential  ellipsoid 
(cf.  Heiskanen  and  Moritz,  1967,  sec.  2-7). 

We  put 


X 

W 


Xo  + aX  • 


U + T 


(3-2) 

(3-3) 


r 
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where  the  differences  6X  = X - X^  and  T = W - U are  considered 

o 

to  be  smallj  T is  called  the  anomalous  potential  (ibid.,  sec. 2- 
13) . 

Thus  (2-5)  becomes 

1 = F(Xq+6X,  U+T)  (3-4) 

and  a Taylor  expansion  gives 

1 = F(Xq,U)  + aT6X  + LT  (3-5) 

plus  higher  order  terms,  which  we  neglect.  Here  a is  the  column 
vector  of  ordinary  partial  derivatives 

*„  ■ !r<xo-u>  <3-«) 

k 

of  F with  respect  to  the  component  X of  the  parameter  vector 
X,  taken  for  the  approximate  values  XQ  and  U;  a is  the  cor- 
responding row  vector,  so  that  aTsX  is  a scalar  product.  The 
term  LT  is  less  elementary:  it  expresses  a linear  operator  L 
acting  on  the  function  T.  The  meaning  of  this  will  be  clear  from 
the  examples  to  follow. 

By  means  of  the  substitution 

61  » F (X , W)  - F(X0,U)  . (3-7) 

the  nonlinear  system  (3-1)  thus  becomes  the  linear  system 

6lx  » a*6X  + LJ  , 

6 1 j - a£6X  + L2T  , 


(3-8) 
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The  linearization  process  will  be  made  clear  by  considering 
some  basic  special  cases. 

Astronomical  and  gravimetric  observations.  The  equations 
to  be  linearized  are  (2-9),  (2-10),  and  (2-11).  Here  we  are  con- 
siderably helped  by  the  fact  that  these  expressions  are  just  or- 
dinary functions  of  Wx,  Wy,  Wz;  X is  simply  the  coordinate 
vector  (2-12). 

Therefore,  we  first  linearize  the  gradient  vector  (2-7). 
Using  index  notation,  we  write  x = y = x2,  z = x3  and 


grad  W 


rw ' 

rw  i 

X 

i 

U 

— 

Wo 

y 

2 

W 

Wo 

z 

3 

(3-9) 


Thus 


W. 


3 W 
3X, 


(i  = 1,2,3)  . 


(3-10) 


X = X, 
O k 


The  derivatives  are  taken  at  the  original  point  with  coordinates 
X = [xj  (k  = 1,2,3)  . (3-11) 

The  approximation  point  is 

(3-12) 

in  obvious  notation, 

xk.x°*sxk.  (3-13) 

Then  (3-10)  becomes 

\ IJ  ^ ( o 2 LI  1 

6 x.  . (3-14) 


aW 

faw  ’ 

x 

f 32W 

3xi 

l»*i( 

T 

0 

(3xi3XjJ 

j 


using  the  summation  convention  (summation  over  the  repeated  index 
j).  The  notation  (')  indicates  that  the  respective  quantity  is 


■v 
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to  be  taken  at  the  approximation  point  (3-12);  is,  of  course, 

taken  at  the  original  point  X. 

We  now  introduce  W = U + T and  obtain 


aw 

fall  1 

X 

raT  1 

X 

r a2u 

6 X + 

0 3 

f 32T 

3xi  = 

isxij 

0 

[sxiJ 

0 

ax  ax 

« •'  < 

3 X 3 X 

l 1 Jj 

,5X3 


(3-lb) 


The  last  term  is  already  of  second  order  ( T and  ax..  are  first- 
order  quantities)  and  will  be  neglected.  We  further  put 


o 


(3-16) 


faT  ) 


ax 


a2u  1 


axiaxj 


M 


ij  * 


Thus  (3-15)  becomes  finally 


Wi  * Ui  + Ti  + • 


(3-17) 

(3-18) 


(3-19) 


completing  the  linearization  of  the  gravity  vector  (3-9). 

The  straightforward  way  to  linearize  equations  (2-9)  to 
(2-11)  is  to  substitute  (3-19)  into  these  equations  and  to  expand 
the  functions  in  the  usual  way  by  Taylor's  theorem,  considering 
the  fact  that  the  second  and  third  term  on  the  right-hand  side  of 
(3-19)  are  small.  This  is  simple  but  laborious;  more  efficient  is 
an  indirect  procedure. 

We  combine  (2-6)  and  (2-7)  into  the  equation 


Wx  ■ -gcosacosA  , 

W2  ■ -gcosasInA  , (3-20) 

W3  ■ -gsln#  ; 

here  all  quantities  refer  to  the  original  point  X. 


l 


- 
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In  an  analogous  way  we  write 


U1  = - yCOS  <}>COS  A 
U 2 = -ycos<f>sinA 
U-  = -ysin$  . 


(3-21) 


Here  all  quantities  refer  to  the  approximation  point:  y is  nor- 
mal gravity,  and  <p  and  A are  normal  latitude  and  longitude. 
Cf.  (Heiskanen  and  Moritz,  1967,  p.315),  where  these  normal  geo- 
graphical coordinates  have  been  denoted  by  and  a*  . 

We  put 


$ = 4>  + 6 4>  , 
A = A + 6A  , 

g = y + 6g  , 


(3-22) 


substitute  into  (3-20)  and  expand  by  Taylor.  The  result  is  readi 
ly  found  to  be 


(3-23) 


U2  + 0 yCO S <f> 6 A 


where  the  matrix 


’sin<j>cosA  sinA  -cos<(>cosa 

Q = sin^sinA  -cosa  -cos<j>sinA 

- co  s <{)  0 - s 1 n <j) 


(3-24) 


is  obtained  by  differentiation  of  (3-21). 

On  the  other  hand  we  have  (3-19),  which  may  be  written  in 
the  form 


'W  ] fu  1 
l l 

W2  = U2  + gradT  + M6X  . 


(3-25) 


->JLJ 
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I 

\ 

i 

i 

i 


The  matrix  Q is  easily  seen  to  be  orthogonal  (why?);  therefore 
its  inverse  is  simply  the  transpose: 

Q"1  = QT  . (3-26) 

Therefore,  the  comparison  of  (3-23)  and  (3-25)  gives 


y 5 ♦ 

Y CO  S $ 6 A 

«g 


QtM6X  + QTgrad  T 


(3-27) 


which  completes  the  linearization  of  astronomical  latitude  $> 
and  longitude  A and  of  measured  gravity  g. 

It  is  evident  that  (3-27)  is,  indeed,  a linear  function  of 
the  components  fix,  fiy,  fiz  of  the  vector  fiX.  As  regards 
QTgrad  T , it  gives  for  each  difference  54,  5A,  fig,  a linear  ex- 
pression of  the  form 


aT  . aT 
“137  + a23y 


+ 


3T 

a33Z 


(3-28) 


The  operation  expressed  by  the  functional  L consists  in  forming 
the  partial  derivatives  and  taking  a linear  combination  of  them. 
Since  differentiation  is  a linear  operation,  L is  indeed  a lin- 
ear functional . 

Direction  and  distance  measurements.  The  straightforward 
approach  is  to  differentiate  equations  (2-14),  (2-15)  and  (2-16), 
as  outlined  In  (Heiskanen  and  Moritz,  1967,  pp. 220-221).  The  re- 
sult will  be  differential  formulas  of  form  of  eq.  (5-83),  ibid. 
The  actual  work  Is,  however,  quite  cumbersome  though  not  diffi- 
cult. 

Again,  an  indirect  approach  might  be  preferable.  We  put 


(3-29) 


t 


i 

l 


* 


ssinZcosA 

ssinZsinA 

scosZ 


- Y . 
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Then  the  vector  Y so  defined  is  related  to  the  difference  vector 


>b  “ * 
ax  = yQ  - 

by  the  linear  transformation 
Y = RaX  , 

where  R is  the  orthogonal  matrix 


(3-30) 


(3-31) 


-si n$cosA 

-si n$s i nA 

COS4> 

R = 

- s i n A 

COSA 

0 

COS4>COSA 

cos$s i nA 

s i n * 

This  is  clear  because 

u,v,w  can 

be  i 

coordi nates 

in  a local 

system  in 

which 

(3-32) 


tion  of  the  gravity  vector  and  the  axes  u and  v point  north 
and  east,  respecti vely ; the  matrix  R is  formed  by  the  components 
in  the  xyz  system,  of  the  unit  vectors  e‘ , eu,  n corresponding 
to  the  uvw  coordinate  axes  (Heiskanen  and  Moritz,  1967  , p p . 2 1 8 - 
219). 

The  differentiation  of  (3-29)  gives,  in  analogy  to  (3-23) 


a 6 Z 

6Y  = S aSinc6A 


(3-33) 


where  S is  the  orthogonal  matrix 

cosecosa  -slna  slnecosa 

S = cosesina  cosa  sintsina 
-sine  0 cose 


(3-34) 
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Here  we  have  designed  by  a ,c,a  the  "normal"  equivalents  of  the 
observables  A,Z,s,  so  that 

A = a + 6A  , 

Z = c + 6Z  , (3-35) 

5 = a + 6S  . 

The  quantities  a,s,o  can  be  computed  from  (2-14),  (2-15),  ana 
(2-16)  by  using  approximate  coordinates  X and  replacing  <t,A 
by  $,x. 

By  differentiation  of  (3-31),  on  the  other  nand,  we  find 

6 Y = R5AX  + 6 RaX  . (3-36) 


The  combination  of  (3-33)  and  (3-36)  gives,  in  view  of  the  ortho- 
gonality of  the  matrix  S, 


a6Z 

asin^aA 

6S 


STR6AX  + ST6RaX  . 


(3-37) 


The  second  term  on  the  right-hand  side  is  easily  found  in 
an  indirect  way.  The  matrix  R,  by  (3-32),  depends  on  # and  a; 
therefore  6R  will  be  a linear  function  of  6*  and  6A.  The 
term  St6RaX  represents,  therefore,  the  effect  of  6*  and  6A 
on  6 Z and  aA  (there  is,  evidently,  no  effect  on  as  because 
s is  independent  of  the  gravity  field);  this  is  nothing  else  but 
the  well-known  effect  of  the  deflection  of  the  vertical  on  azi- 
muth A and  zenith  distance  Z. 

The  effect  on  the  zenith  distance  Z is 


3 Z =»  tcosa  + nsina 


(3-38) 


(Heiskanen  and  Moritz,  1967,  p.190,  eq .( 5-20))  and  on  the  azimuth 


_i — l 
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3 A = 5 s i nacot;;  + n ( t a n <{>  - cosacotc)  (3-39) 

(ibid.,  p.186,  eq. (5-13) ) . We  have  used  the  symbols  aZ  and  sA 
to  indicate  the  partial  influence,  on  Z and  A,  of  the  changes 
6<t>  and  6 A , which  are  related  to  the  deflection  components  j; 
and  n by 

£ = 64>  , n = 6 A CO  S <f>  . (3-40) 

The  comparison  of  (3-24)  and  (3-32)  shows  that,  for  <j>  = 4, 
and  A = x, 

R = -QT  i (3-41) 

the  geometrical  interpretation  of  this  fact  is  left  to  the  reader. 

In  view  of  the  relations  (3-38)  to  (3-41),  eq.(3-37)  takes 
the  final  form 


o 6 Z 

X - 6 X ' 

osinsSA 

= -stqt 

Q P 

% - % 

+ kT6$ 

[_C  0 S 6 A_ 

. «s  . 

6 Z - 6 Z_ 

L Q PJ 

where 


K 


a C0Sa 
os i nacoss 
0 


osina 

a(tan<(.sin£  - cosacos s) 
0 


(3-42) 


(3-43) 


the  matrices  Q and  S are  given  by  (3-24)  and  (3-34). 

These  examples  will  illustrate  how  the  linearized  equations 
(3-8)  can  be  obtained.  Similar  linearizations  can  be  found  in 
(Eeg  and  Krarup,  1975)  and  (Grafarend,  1977). 

Satellite  observations.  Satellite  observations  can  be  lin- 
earized in  the  same  way,  as  outlined  In  (Kaula,  1967, p. 67)  or 
(Heiskanen  and  Moritz,  1967,  pp.352  - 355).  This  Is  theoretically 
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straightforward  but  computationally  very  laborious.  For  direction, 
distance,  and  Doppler  measurements,  linearizations  can  be  found 
in  many  forms  in  the  literature.  Therefore,  we  need  not  go  into 
details  here.  We  only  mention,  that  the  1 i near  functional  LT  is, 
for  satellite  observations,  usually  expressed  in  spherical  har- 
monics: 


LT 


? ? (a  6C  + b 6S  ) 

£,_  nm  nm  nm  nm 
n*  2 m=0 


(3-44) 


which  is  a linear  combination  of  the  spherical -harmonic  coeffi- 
cients 6 C and  6 S of  T. 
nm  nm 

For  second-order  gradients  the  linearization  is  straight- 
forward: we  have 


3 2 V 

r 3 2 v l 

E 

3 3 V 

E 

BX^X 

lsvx3J 

0 

3 X 3 X 3 X 

1 j 

which  is  of  form 

(3-5). 

VE 

being  the 

«X  + 

o k 3xi3Xj 


(3-45) 


vitational  potential.  The  vector  6xk  can  be  further  expressed 
in  terms  of  the  orbital  parameters;  cf.  eq.(9-35)  of  (Heiskanen 
and  Moritz,  1967,  p.352). 

A linearization  for  satel 1 i te-to-satel 1 i te  tracking  can  be 
found  in  (Kryfiski,  1978). 


4.  Determination  of  the  Gravity  Field 
as  an  Improperly  Posed  Problem 

A problem  is  called  properly  posed  if  the  solution  satis- 
fies the  following  three  requirements: 

(1)  existence, 

(2)  uniqueness, 

(3)  stability. 

This  means  that  a solution  must  exist  for  arbitrary  (within  a 
certain  range)  data,  that  there  must  be  only  one  solution,  and 
that  this  solution  must  depend  continuously  on  the  data.  If  one 
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or  more  of  these  requirements  are  violated,  then  we  have  an 
improperly  posed,  or  i 1 1 - p o s e d , probl em. 

For  a long  time  it  was  thought  that  only  properly  posed 
problems  are  physically  meaningful.  In  fact,  deterministic  pro- 
cesses, as  considered  in  classical  mechanics,  depend  uniquely 
and  continuously  on  the  initial  data--this  is  the  essence  of 
causal i ty--and  thus  correspond  to  properly  posed  problems. 

Only  relatively  recently  is  was  recognized  that  there  are 
important  problems  that  are  not  properly  posed.  There  is  now  an 
extensive  literature  on  improperly  posed  problems;  we  mention 
only  two  easily  accessible  books:  (Lavrentiev,  1967)  and,  espe- 
cially, (Tikhonov  and  Arsenin,  1977),  and  the  review  article 
(Nashed,  1974).  Geodetic  applications  are  considered  in  (Schwarz, 
1978);  for  instance,  the  downward  continuation  of  gravity  is  an 
ill-posed  problem.  Also  the  work  by  Neyman  (1977)  should  be  men- 
tioned. 

Our  present  task,  the  determination  of  the  earth's  gravi- 
tational field  from  measurements,  is  a typical  improperly  posed 
problem.  The  potential  is  so  irregular  that  it  cannot  be  comple- 
tely described  by  any  finite  set  of  parameters;  on  the  other  hand, 
we  have  only  a finite  number  of  measurements.  Hence,  there  is  no 
unique  solution,  and  Condition  2 is  violated. 

Using  the  standard  notation  of  (Tikhonov  and  Arsenin,  1977), 
we  may  write  equations  (3-8)  in  the  form 

Az  = u , (4-1) 


if  we  put 
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A is  a linear  operator,  expressing  the  fact  that  equations  (3-8) 
permit  the  computation  of  the  observations  6lk  if  <5X  and  T 
are  supposed  to  be  given;  the  operations  involved  in  these  com- 
putations are  clearly  linear.  The  operator  A,  therefore,  com- 
prises the  vectors  a£  and  the  linear  functionals  Lk. 

In  mathematical  terms,  6X  is  a vector  belonging  to  m- 
dimensional  Euclidean  space  R , and  T can  be  considered  a mem- 
ber of  a Hilbert  space  H of  harmonic  functions;  u belongs  to 
R . Therefore,  the  linear  operator  A (4-1)  defines  a linear 
mapping: 


A:  R x H -►  R . (4-4) 

m q ' 

The  solution,  if  it  exists,  may  be  written  formally  in  the  form 

z = A’1 u , (4-5) 

but  it  will  certainly  not  be  unique.  Therefore,  A-1  is  not  an 
inverse  operator  in  the  usual  sense;  it  has  the  character  of  a 

generalized  inverse  operator  (analogous  to  generalized  matrix  in- 

.1 

verses). 

For  the  solution  of  our  problem  we  shall  try  to  apply  stan- 
dard mathematical  techniques  for  Ill-posed  problems.  Nashed  (1974, 
p.295)  mentions  the  following  possible  approaches: 

(a)  a change  of  the  concept  of  a solution; 

(b)  a change  of  the  spaces  and/or  topologies; 

(c)  a change  of  the  operator  Itself; 

(d)  the  concept  of  regul arization  operators; 

(e)  probabilistic  methods  or  well-posed  stochastic  exten- 
sions of  111-posed  problems. 

These  approaches  may  overlap  In  various  ways. 

We  shall  especially  use  the  approaches  (b),  (d),  and  (e); 
they  will  be  considered  In  the  following  sections  5,6,  and  7. 

We  finally  point  out  that  the  vector  X will  only  have 
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linearly  independent  components  and  that  also  all  equations  of 
the  system  (3-8)  are  supposed  to  be  linearly  independent.  In 
other  terms,  we  shall  work  with  systems  (4-1)  that  have  full  rank. 


5 . Pure  Collocation 

Let  us  suppose  that  the  systemat i c parameters  X (coordi- 
nates, etc.)  are  known  with  sufficient  accuracy.  We  then  have 
6X  = 0,  and  the  system  (3-8)  reduces  to 

V * J 1,2 q ; (5’1} 

for  simplification,  we  have  replaced  6lR  by  1^.  In  order  to 
find  a solution,  we  apply  Approach  (b)  mentioned  at  the  end  of 
the  preceding  section:  a change  of  the  solution  space. 

We  shall  try  to  approximate  the  desired  function  T by  a 
linear  combination  f of  suitable  linearly  independent  base 
functions  <t>.  , <t>„ , . . . , 4>  : 

12  q 

T(P)  ± f(P)  =•  ? b ♦ (P)  , (5-2) 

k=l  K 

P denoting  the  space  point  at  which  these  functions  are  being 
considered,  and  bR  denoting  suitable  coefficients.  Since  T is 
harmonic  outside  the  earth's  surface,  the  base  functions  4>k 
must  also  be  harmonic  functions. 

Since  there  are  q independent  functions  $k,  the  space 
of  linear  combinations  (5-2)  is  q -dimensional,  so  that  the  so- 
lution space,  as  well  as  the  observation  space,  is  q -dimension- 
al, and  the  operator  A,  consisting  of  the  q functionals  Lk> 
must,  for  this  particular  case,  reduce  to  an  q x q matrix,  which 
in  general  has  a regular  Inverse;  hence  the  problem  will,  in 
general,  become  well-posed. 

Note  that  this  Is  made  possible  by  changing  the  solution 
space  from  Inflnltely-dlmenslonal  Hilbert  space  to  the  q -dimen- 
sional space  of  linear  combinations  (5-2). 


j~-  ■’  - 
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For  the  special  case  of  interpolation,  details  can  be  found 
in  (Moritz,  1978b,  c). 

We  substitute  (5-2)  into  (5-1)  and  get,  because  of  lineari- 
ty, 


?■> 


kLj*k 


The  quantities 


(5-3) 


L 


(5-4) 


being  the  values  of  linear  functionals,  are  constants.  Hence  we 
get  the  system  of  q linear  equations 


f A 

k-1 


jkbk 


* 


(5-5) 


which  can  be  solved  for  the  q coefficients  bR  provided  the 
determinant  of  the  matrix  A,.  is  non-zero.  The  substitution  of 
the  bk  into  (5-2)  gives  the  desired  solution. 

The  approximation  of  a function  by  fitting  a linear  com- 
bination of  base  functions  to  a number  of  linear  functionals  is 
col  1 ocation  in  the  sense  of  approximation  thdory,  cf.(Collatz, 
1966,  pp.29). 

As  a matter  of  fact,  interpolation  Is  a special  case  of 
collocation,  when 


V • Kfj)  . (6-6) 

that  is,  when  L ^ associates  to  a function  its  value  at  a parti- 
cular point  P ^ ; this  is  the  so-called  evaluation  functional. 

Kernel  functions.  Consider  now  a positive-definite  symme- 
tric function  K(P,Q) , harmonic  as  a function  of  both  P and  Q. 
The  positive  definiteness  of  this  function  means  that  for  any  N, 
the  N x N matrix  of  elements  K(P^,Pk),  for  any  points  PJt  P2, 
....  PN,  is  positive  definite.  Harmonicity  means  that  Laplace's 


T 
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equation  is  satisfied  both  at  the  point  P (with  Q held  fixed) 
and  at  Q (for  f i xed  P) : 

(5-7) 
(5-8) 

the  domain  of  harmonicity  may  be  the  exterior  of  a certain  sphere 
completely  inside  the  earth. 

Such  a function  K(P,Q)  will  be  called  a kernel  function. 
Assume  now  that  base  functions  <t> ^ have  the  special  form 

*k(P)  = LkK(P’Q)  * (5-9) 

where  L®  means  that  the  linear  functional  is  applied  to  the 
variable  Q. 

Then  the  matrix  (5-4)  takes  the  form 

cjk  = L^L2(K,Q)  , ; (5-10) 


ApK(P,Q)  = 0 
AqK(P.Q)  = 0 
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and  (5-5)  becomes 

j.Vi.-V  <5-,1> 

We  solve  this  linear  system  and  substitute  the  result  into  (5-2), 
which  by  (5-9)  has  the  form 


f(P)  * ? b ♦ (P)  = ? b L®K(P ,Q)  . (5-12) 

j-l  J 3 j-l  3 J 

The  result  may  be  written  as 
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using  the  abbreviation 


Cpj  = L*K(P,Q)  . 


(5-14) 


This  is  the  basic  formula  for  analytical  collocation  with 
kernel  functions. 

Interpolation  with  kernel  functions  is  the  special  case 
(5-6),  so  that  simply 


*j(P)  - K(P,P.)  = Cpj  . 


■ K(Pj*PJc) 


(5-15) 

(5-16) 


This  interpolation  method  is  known  from  mathematics  (Meschkowski , 
1962,  p . 114) i collocation  with  kernel  functions  has  been  applied 
in  geodesy  first  by  Krarup  (1968,  1969). 

A remark  on  terminology:  the  approximating  function  is  ob- 
tained by  fitting  the  given  functionals  exactly  (measuring  errors 
are  not  taken  into  account);  therefore  one  speaks  of  "exact'1  or 
"pure"  collocation.  Sometimes  it  is  also  called  "deterministic", 
to  distinguish  It  from  stochastic  methods  (Dermanls,  1976,  p.56; 
Moritz,  1976,  p.25),  but  the  naiiie  "exact",  or  "pure",  or  "analy- 
tical", collocation  seems  to  be  better  because  the  concept  "deter- 
minism" already  has  a well-established  meaning,  namely  causality 
in  the  sense  of  classical  physics. 

Minimum  norm  property.  The  solution  (5-13)  has  an  impor- 
tant property:  among  all  possible  solutions  (In  some  Hilbert  spa- 
ce of  harmonic  functions)  of  the  system  (5-1),  the  solution  (5-13) 
has  the  minimum  norm: 


||  f||  * minimum  , 


(5-17) 


If  the  norm  Is  defined  by  the  Inner  product 

II  fll  2 “ (f.f) 


(5-18) 
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in  the  Hilbert  space  with  the  kernel  function  K(P,Q),  so  that 

there  holds  the  reproducing  kernel  property  > I 

(f(P).K(P.Q))  - f (Q)  , (5-19)  | 

where  (•,«)  denotes  the  inner  product  with  respect  to  P (Mesch- 

kowski,  1962,  p.115;  Krarup,  1969,  p.39;  Tscherning,  1975;  Moritz,  ‘ 

1978b,  p.36). 

1 

6 . Application  of  Tichonov  Regularization 

Since  we  shall  use  the  letter  A for  a different  purpose 
later  in  this  section,  we  shall  denote  the  operator  in  (4-1)  by 
G and  write 

Gz  = u . (6-1) 

Tichonov's  regularization  method  consists  in  minimizing 
the  nonlinear  functional 

Ma[z,u]  = ||  Gz  - u ||  2 + aft ( Z ) , (6-2) 

where  a is  a numerical  parameter  and  n(z)  is  a so-called 
stabilizing  functional  (Tichonov  and  Arsenin,  1977  , pp. 51,57), 
which  may  be  taken  as  the  square  of  some  norm. 

fi(z)  » II  z||  2 (6-3) 

(ibid. , p.72).  In  this  way,  a unique  solution  can  usually  be  ob 
tai ned . 

Using  (6-3)  and  very  slightly  generalizing  (6-2)  by  the 
introduction  of  a second  numerical  parameter  8,  we  get  the  con 
di tion 
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T 


“II  zl|2  + ell  Gz 


u 1 1 * minimum 


(6-4) 


Obviously,  Gz  - u is  the  error  in  satisfying  (6-1),  so 
that  ||  Gz  - u 1 1 will  be  the  error  norm.  Therefore,  (6-4)  mini- 
mizes a linear  combination  of  the  norm  of  the  solution  z and 
of  the  error  norm.  Depending  on  the  relation  between  a and  e, 
a stronger  weight  is  given  to  one  or  the  other  of  these  two  norms. 

The  linear  operator  equation  (6-1)  is,  in  our  case,  nothing 
else  than  the  system  (3-8).  Let  us  simplify  the  notation  by  re- 
pl aci ng 
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61k  by  \ • 


Then  (3-8)  becomes 


\ = a?  + LiT 
12  - a*X  ♦ L2T 


1 = aAX  + L T 

q q q 


6 X by  X 


(6-5) 


(6-6) 


where  X is  a m -vector  (an  m x 1 matrix).  We  finally  put 
’1 


1 


1. 


1 

q 


aT 

al 

aT 

a2 


aT 
L <IJ 


(6-7) 


and 


ll*j 


(6-8) 


Here  1 is  a q -vector  (a  q x 1 matrix),  A is  a q x m ma- 
trix, and  B Is  a linear  operator,  formed  of  the  q linear 

functionals  L,  . 

k 
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With  these  notation,  (6-6)  becomes 


1 = AX  + BT 
which  is  (6-1)  with 


(6-9) 


u 

z 

G 


(6-10) 

(6-11) 

(6-12) 


we  have  used  for  linear  operators  the  same  mode  of  partitioning 
as  for  matrices.  In  fact,  (6-1)  becomes  on  partitioning 


[A  B] 


X 

T 


= 1 


(6-13) 


which  is  (6-9) . 

What  is  Gz  - u in  our  present  case?  We  have 


AX  + BT  - 1 


(6-14) 


which  is  the  amount  by  which  the  exact  equation  (6-1)  or  (6-9) 
is  not  satisfied.  Let  us  denote  this  "(Disclosure"  by  -n: 

AX  + BT  - 1 = -n  . (6-15) 

On  rearrangement  this  becomes 

1 = AX  + BT  + n ; (6-16) 

here  n may  be  interpreted  as  the  effect  of  measuring  errors; 
we  call  n the  "noise".  It  should  be  repeated  that  this  equa- 
tion is  nothing  but  the  system  (3-8),  obtained  by  linearizing 
the  nonlinear  functional  equations  (3-1). 


I 

\ 

» 

I 

l 

f 

r 


* 


32 


Then,  in  (6-4), 

||Gz  - U||  = ||  -n||  = ||  nj|  , 

which  is  the  norm  of  the  q -vector  n.  Any  regular  quadratic 
norm  in  q -dimensional  vector  space  can  be  written 

||  n ||  2 = nTQn  (6-13) 

with  a positive  definite  regular  symmetrix  q x q “weight  matrix" 

Q . Let  us  denote  its  inverse  by  D,  then  (6-18)  becomes 

||  n ||  2 = nTD_1  n . (6-19) 

If  n are  random  quantities  in  a statistical  sense,  then 
D may  be  considered  as  the  covariance  matrix  of  the  measuring 
errors  n.  This  statistical  Interpretation  of  n is  not  neces- 
sary--we  may  consider  the  norm  (6-18)  as  a metric  In  a purely 
geometric  sense--but  it  gives  a clue  as  to  the  proper  choice  of 
the  metric  for  the  vector  X.  If  the  variance  of  a random  quan- 
tity is  small,  then  this  quantity  can  vary  only  within  narrow  li- 
mits. The  larger  the  variance,  the  larger  variations  are  possible; 
and  if  the  variance  goes  to  infinity,  the  variation  of  our  quan- 
tity becomes  completely  free. 

Since  we  allow  our  parameter  X to  vary  freely  and  inde- 
pendently, each  component  should  have  an  Infinite  variance,  or 
zero  weight.  Therefore,  In  an  expression  such  as  (6-18),  with  X 
Instead  of  n,  there  should  be  Q * 0,  which  gives 

||  X||  - 0 . (6-20) 

This  will  be  our  choice  for  the  norm  of  X. 

Finally,  the  norm  for  T will  be  selected  a norm  in  a 
Hilbert  space  with  kernel  function  K(P,Q),  that  Is,  defined  by 
(5-18)  and  (5-19); 
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||  T 1 1 2 = (T»T)  . (6-21) 

Then,  assuming  X and  T independent  of  each  other: 

II  z II  2 - II  X||2  + ||  T 1 1 2 = II  T II  2 . (6-22) 

Therefore,  the  Tichonov  condition  (6-4)  becomes 

oc  ( { T 1 1 2 + 0)|  n 1 1 2 = minimum  (6-23) 


or 


a ( T , T ) + BnTD  *n  = minimum  . (6-24) 

Pure  collocation.  As  a preparation,  let  us  assume  error- 
less observations.  Then  n = 0 and  (6-24)  reduces  to 

(T,T)  = minimum  (6-25) 

(we  have  put  a = 1 without  loss  of  generality).  We  furthermore 
assume  X = 0 (no  systematic  effects).  Then  (6-16)  gives 

BT  = 1 (6-26) 

where  1 is  given.  The  desired  T is  that  function  T,  satis- 
fying (6-26),  which  minimizes  (6-25). 

We  solve  the  problem  by  means  of  a Lagrange  multiplyer. 
Instead  of  minimizing  (6-25)  under  the  side  condition  (6-26),  we 
form  the  unconditional  minimum  of  the  function 

* = — (T  ,T ) - kT(BT  - 1)  , (6-27) 

2 

where  the  q -vector  k serves  as  a Lagrange  multiplyer. 

A necessary  condition  for  a minimum  is  the  vanishing  of 
the  differential  of  »: 
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d$  = (T,dT)  - kTBdT  = 0 . (6-28) 

||  (We  can  form  this  differential  as  if  T were  a vector. 

To  avoid  misunderstanding,  we  point  out,  however,  that  dT  is 
not  an  ordinary  differential  of  T,  but  what  is  called,  in  the 
calculus  of  variations,  a first  variation,  that  is,  a change  in 
the  function  T.  In  fact,  (6-24)  is  a variational  principle,  and 
(6-28)  is  the  corresponding  Euler  equation  for  our  special  case. 
Eq.  (6-28)  may  be  found  by  replacing,  in  (6-27),  T by  T + ex, 
where  e is  a small  parameter: 

= -^-(T  + ET »T  + ex)  - kT(BT  + eBx  - 1)  = 

- ?(T,T)  - kT(BT  - 1)  + Je(T,x)  + £e(x,T)  - 

- ekTBx  + ^ e2(  x,  x)  . 

By  symmetry,  (t,T)  = (T,x).  We  subtract  (6-27)  and  divide  by  e. 
On  letting  e -*■  0,  we  thus  get 

lim = (T,t)  . kTBx  „ o , (6-29) 

e-*o  *■ 

which  is  (6-28),  with  dT  = ex.) 

The  function  dT  in  (6-28)  is  completely  arbitrary;  it 
need  not  even  be  small  since  a numerical  factor  does  not  matter. 
(Of  course,  dT  must  belong  to  the  Hilbert  space  under  consider- 
ation. ) 

By  the  reproducing  property  (5-19), 

dT(Q)  = (dT(P),K(P,Q)) 
or,  briefly, 

dT  » (dT , K)  « ( K,dT ) . 


(6-30) 

(6-31) 


Hence , 


BdT  = (BK.dT) 
and  (6-28)  becomes 


(6-32) 


(T,dT)  - ( kTBK , dT ) = 0 


or 


(T  - kTBK,dT)  = 0 . (6-33) 

Since  dT  is  arbitrary,  there  must  be 
T - kTBK  = 0 


T = kTBK  . (6-34) 

This  is  an  important  result.  What  does  it  mean?  In  view 
of  (6-8)  and  (6-30),  this  is  nothing  else  than 

T(Q)  = ! k,L^K(P,Q)  i (6-35) 

j = l 3 -> 

means  the  operator  l-r  applied  to  the  variable  P.  Now,  (6- 
35)  is  identical  to  (5-12),  with  P and  Q interchanged  and 
b = k...  Thus,  the  best  approximation  for  T(Q)  j_s,  in  fact,  a 
linear  combination  of  the  base  functions  (5-9) ! 

The  rest  is  straightforward.  Considering  T a scalar,  we 
may  transpose  (6-34): 

T = (BK)Tk  , (6-36) 

and  substitute  into  (6-26): 

B(BK)Tk  = 1 . 


(6-37) 
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The  q x q matrix 

C = B(BK)t  (6-38) 

has,  by  (6-8),  the  elements 

Cij  = L*l2K(P,Q)  . (6-39) 

Now  (6-37)  may  be  solved  for  k: 

k = C_11  , (6-40) 

so  that  (6-36)  becomes 

T = ( BK)tC“ 1 1 . (6-41) 


In  view  of  (5-14),  this  is  identical  to  (5-13). 

We  thus  have  obtained  (5-13)  as  a consequence  of  the  mini- 
mum norm  principle  (6-25).  There  are  shorter  and  more  complete 
proofs  (the  condition  (6-28)  is  necessary  but  not  sufficient); 
the  advantage  of  the  present  derivation  is  the  treatment  as  a 
straightforward  solution  of  a variational  principle  by  standard 
techniques  (Euler  equation).  Furthermore,  it  will  essentially 
simplify  .he  treatment  of  the  general  case. 

The  case  a * e ■ 1 . As  a second  step,  let  us  consider 
the  general  equation  (6-16),  but  put,  in  the  Tichonov  condition 
(6-24),  o * B » 1,  so  that 

(T,T)  + nTD_1n  « minimum  , (6-42) 

to  be  solved  under  the  side  condition  (6-16).  We  thus  have  to 
form  the  unconditional  minimum  of  the  function 

* - ?(T,T)  + ^nTD*1n  - kT (AX  * BT  + n - 1 ) . (6-43) 

| 
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The  differential  is 

do  = (T,dT)  + nTD_1dn  - kT(AdX  + BdT  + dn)  , (6-44) 

where  dX  and  dn  are  ordinary  vector  differentials.  On  re- 
arranging we  get,  using  (6-32), 

d<t>  = (T  - kTBK , dT ) + ( nTD~ 1 - kT)dn  - kTAdX  = 0 . (6-45) 

Since  dT,  dn,  and  dX  are  arbitrary,  do  = 0 can  only  hold  if 


T - k1 BK  = 0 , 

nTD~ 1 - kT  = 0 , 

kTA  = 0 . 


The  first  equation  gives 


(6-46) 

(6-47) 

(6-48) 


T = k BK  , 


(6-49) 


identical  to  (6-34)  or  (6-35).  Again,  the  solution  is  a linear 
combination  of  base  functions  (5-9)! 

The  transpos i tion  of  (6-49)  gives 


T = (BK)k  , 


(6-50) 


so  that 


BT  = B(BK)Tk  = Ck  , 


using  the  abbreviation  (6-38).  Eq.(6-47)  gives 


(6-51) 


nT  = kTD 


or 


n = Ok  . 


(6-52) 


Eq.(6-16)  may  be  written 


38 


1 - AX  = BT  + n , (6-53) 

and  substituting  (6-51)  and  (6-52)  we  get 

1 - AX  = (C  + D)k  , (6-54) 

so  that 

k = (C  + D ) — 1 ( 1 - AX).  (6-55) 

We  substitute  this  into  (6-48),  transposed  as 
ATk  = 0 , 
obtaining 

AT(C  + D ) ~ 1 1 - At(C  + D ) " 1 AX  = 0 , 

so  that 

X * [at(C  + D ) ” 1 a]  1 At(C  + D)_11  , (6-56) 

which  determines  the  parameter  vector  X.  The  substitution  of 
(6-55)  into  (6-50)  then  gives  the  potential: 

T - (BK)t(C  + D ) ” 1 ( 1 - AX)  . (6-57) 

The  general  case.  Take  finally  the  general  Ttchonov  condi- 
tion (6-24) 

a (T ,T ) + 0nTD-1n  * minimum  (6-58) 

for  solving  the  equation  (6-16) 


AX  + BT  + n - 1 . 


(6-59) 
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The  condition  (6-58)  is,  for  a t 0,  equivalent  to 

i 

(T,T)  + nT(|-D)-1n  = minimum  , (6-60) 

so  that  we  only  have  to  replace,  in  (6-56)  and  (6-57),  D by 
aD/6.  This  gives 

X = [AT(BC  + a D ) - 1 a]  1AT(6C  + aD ) “ 1 1 , (6-61) 

T = (BBK)T(BC  + aD ) ~ 1 ( 1 - AX)  , (6-62) 

which  is  the  solution  of  (6-59)  under  the  general  Tichonov  con- 
dition (6-58)  . 

Obviously,  to  various  ratios  a : 3 there  corresponds  a 
different  weighting  between  the  square  of  the  "function  norm", 

(T,T),  and  of  the  "error  norm",  nTD_1n.  Pure  collocation,  with 
the  condition  (6-25)  and  n = 0,  fits  the  solution  exactly  to 
the  data.  This  is  unsuited  for  real  data  in  the  presence  of  mea- 
suring errors,  because  then  the  solution  is  distorted  by  faith- 
fully reproducing  all  measuring  errors;  we  risk  to  get  spurious 
oscillations  (Eeg  and  Krarup,  1975,  p . 1 1 1 ) ^ ' 

Therefore,  some  balance  between  a and  3 in  (6-58)  must 
be  found.  But  how?  We  might  use  some  tri al -and-error  procedure, 
but  this  does  not  seem  very  satisfactory.  A theoretically  moti- 
vated, in  a certain  sense  optimal,  solution  to  this  problem  is 
found  by  statistical  considerations,  as  we  shall  see  in  the  fol- 
lowing section. 

In  their  "Integral  Geodesy",  Krarup  and  Eeg  (1975)  give 
equations  equivalent  to  (6-61)  and  (6-62),  but  with  different 
weighting  in  the  two  equations;  In  (6-61)  they  use  6 = 1 - a 
(ibid. , p . 119,  their  eq . ( 5 ) ) and  in  (6-62)  they  use  3 ■ a = 1 
(ibid.,  p . 1 12 , their  eq.(14)).  This  follows  from  employing  two 
different  minimum  principles.  It  seems,  however,  preferable  to 
derive  both  equations  from  the  same  variational  principle  (6-58). 

We  also  mention  that  the  result  for  the  "errorless"  con- 
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dition  (6-25)  cannot  be  obtained  by  simply  putting  a = 1,  e * 0, 
as  might  be  expected  at  first  sight.  The  essential  feature  with 
(6-25)  is  that  n = 0 and  D = 0.  We,  therefore,  have  to  put 
D = 0 in  (6-61)  and  (6-62).  As  a consequence,  0 cancels  then, 
and  we  obtain 

X = (AtC"1A)'1AtC"11  , (6-63) 

T = ( BK)tC_ 1 ( 1 - AK)  . (6-64) 

For  A = 0 (no  systematic  parameters),  the  last  equation  reduces 
to  (6-41),  as  it  should. 

A final  word  on  the  solution  of  these  variational  princi- 
ples by  an  Euler  equation.  Any  Euler  equation  gives  only  a neces- 
sary, not  a sufficient,  condition  for  a minimum.  It  is  not  too 
difficult  to  show  that  our  solutions  do  indeed  give  a minimum. 

The  proof  goes  along  well-known  lines  (cf.  Moritz,  1972,  pp. 22-23). 


7 . Least-Squares  Collocation 

One  of  the  possible  approaches  to  Improperly  posed  problems 
mentioned  at  the  end  of  sec. 4 are  stochastic  methods.  We  have  al- 
ready been  close  to  such  an  approach  in  the  preceding  section, 
when  we  Interpreted  the  matrix  D In  (6-19)  as  the  covariance 
matrix  of  the  measuring  errors;  however,  this  was  not  essential 
since,  from  an  analytical  point  of  view,  we  could  as  well  have 
used  any  other  positive-definite  regular  square  matrix. 

A more  thoroughgoing  stochastic  approach  consists  In  the 
attempt  to  Interpret  the  norm  ||  T||  statistically,  as  well  as 
the  error  norm  ||  n||  . This  can  be  achieved  by  formally  consider- 
ing the  anomalous  gravity  potential  T as  a stochastic  process. 
This  could  be  a purely  formal  "stochastic  extension",  motivated 
by  the  Irregularity  of  the  anomalous  gravity  field. 
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Historically,  this  was  the  way  that  led  to  collocation  in 
its  geodetic  application.  The  statistical  treatment  of  the  gra- 
vity anomalies  was  initiated  by  DeGraff -Hunter  in  1935,  but  did 
not  find  much  response.  The  contemporary  rapid  development  start- 
ed from  Hirvonen's  work,  published  in  1956:  the  explicit  treat- 
ment of  the  anomalous  gravity  field  as  a stochastic  process  in 
1959  by  Kaula,  the  geodetic  application  of  least-squares  predic- 
tion techniques  of  stochastic  processes  by  Moritz  in  1962,  genera- 
lized by  Kaula  in  1963,  and  the  extension  to  the  estimation  of  T 
from  given  functionals,  which  is  collocation  in  the  proper  sense, 
by  Krarup  in  1968-9. 

Stochastic  process  techniques  form  a very  convenient  tool 
and  an  intuitive  terminology  (variances,  covariance  functions, 
etc.),  and  they  permit  the  estimation  of  accuracies,  which  is  of 
particular  importance  for  geodetic  applications. 

There  is  no  fundamental  objection  against  embedding  the 
actual  Earth,  with  its  potential  T,  in  an  ensemble  of  "sample 
earths",  making  possible  the  application  of  stochastic-process 
techniques.  Similar  approaches  are  frequently  applied  (conscious- 
ly or  unconsciously)  in  various  other  fields. 

From  the  point  of  view  of  logical  simplicity,  however,  it  t- 

appears  preferable  to  find  an  interpretation,  which  retains  the 
convenient  formal  apparatus  of  stochastic  processes,  but  is  re- 
stricted to  our  Earth  only,  without  introducing  fictitious  other 
"earths".  t 

This  is  made  possible  by  an  application  of  Norbert  Wiener's 
"covariance  analysis  of  individual  functions",  in  our  case,  of 
the  actual  anomalous  gravity  potential  T (Moritz,  1972, sec. 8). 

In  this  case,  the  covariance  function  of  T is  defined  as 

C(P,Q)  - M1{TpT0>  , (7-1) 

as  a homogeneous  and  isotropic  average  M1  of  the  product  tpTq» 
where 

y 

. . ... 
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M { •}  = A /2,r  /*  /2n(  • )sinededxda  (7-2) 

1 01T  x =o  e=o  a=o 

is  an  average  over  all  points  (e,x)  and  all  azimuths  a (ibid., 
p . 1 1 4 ) , which  can  be  considered  as  an  average  over  rotation  group 
space , 

If  we  consider  rotation  group  space  as  a probability  space 
and  on  this  space  a stochastic  process,  the  sample  functions  of 
which  differ  only  by  rotations,  then  the  probability  average  co- 
incides,  by  definition,  with  the  "space  average"  (7-1),  so  that 
our  process  will  be  (trivially)  ergodi c (Moritz,  1972,  p.119). 

This  is  elaborated  in  detail  in  (Moritz,  1978a),  where  ergodic 
processes  on  the  sphere  are  discussed,  They  are  of  necessity  non- 
Gaussian,  since  Lauritzen  has  shown  in  1971  that  ergodic  Gaussian 
processes  on  the  sphere  are  impossible. 

The  statistical  expectation  in  the  probability  space  of 
the  measuring  errors  n are  denoted  by  M2> 

Thus  the  measuring  errors  may  be  considered  as  "physically 
stochastic"  quantities.  The  anomalous  potential  T is  only  for- 
mally a space  average.  If  we  understand  statistics  as  the  study 
of  a large  amount  of  data  and  their  average  properties  (for  in- 
stance, an  analysis  of  the  global  human  population),  then  we  may 
say  that  our  treatment  of  T Is  statistical , though  T itself 
is  not  a stochastic  quantity  In  a strictly  physical  sense  (Mo- 
ritz and  Sanso,  1978). 

Least-squares  collocation  is  defined  as  that  collocation 
method  which  minimizes  the  variance  of  the  estimated  quantities, 

In  our  case,  of  the  potential  T and  the  parameter  vector  X. 

Let  T*  be  the  true  potential  and  T its  estimate  from  the  gi- 
ven data,  affected  by  measuring  errors,  at  a particular  point  P. 
Then  the  variance  of  the  estimate  T , or  the  square  of  the  pre- 
diction error  in  , is  defined  as 

p 

"p  * M1M2{<TP  * V2>  • (7“3) 
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that  is,  as  the  average  of  the  square  of  the  true  prediction  er- 
ror  T - T at  some  point  P,  extended  over  all  measuring  er- 


rors 


(M2) 


and  all  points  of  the  sphere  (M  ) . 

This  definition  is  very  natural;  it  corresponds  to  a glo- 
bal average. 


Similarly  the  variance  m^ 
nent  X.  of  the  parameter  vector 


mp  = minimum  , 


nr  = minimum  » 


of  the  estimate  for  a compo- 
X is  defined.  The  conditions 

(7-4) 

(7-5) 


for  all  points  P of  the  sphere  and  all  components  of  X,  define 
least-squares  collocation. 

Special  cases  are  least-squares  prediction  of  gravity, 
which  uses  only  (7-4)--cf.  (Heiskanen  and  Moritz,  1967  , p . 268) -- 
and,  of  course,  least-squares  adjustment,  which  uses  only  (7-5). 

The  result  of  conditions  (7-4)  and  (7-5)  are  as  follows: 
we  get  collocation  in  a Hilbert  space,  in  which  the  kernel  func- 
tion is  given  by  the  covariance  function  (7-1),  and  T and  X 
are  then  expressed  by  (6-56)  and  (6-57).  Putting 


C + D = C = Z , 

XX 

BK  = C 

sx 


(7-6) 

(7-7) 


we  get  the  well-known  formulas  of  least-squares  collocation 


= (AaTT”1A)-1At£’"1 
-l 


= c Z 

sx 


(1  - AX) 


(7-8) 

(7-9) 


Here  we  have  replaced  T by  s,  because  the  anomalous  potential, 
or  any  other  quantity  of  the  anomalous  gravitational  field,  is 
also  called  a "signal",  denoted  by  s.  A pleasant  consequence  of 
least-squares  methods  with  respect  to  linear  transformation  is 
that  an  equation  of  form  (7-9)  holds  for  the  direct  estimation 


4 m-  • 
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of  any  anomalous  field  quantity,  all  these  quantities  (signals) 
being  linear  functions  of  T.  In  fact,  eq.(2-36)  of  (Moritz, 
1972),  has  been  derived  for  the  estimation  of  any  signal. 

Thus,  least-squares  collocation  corresponds  to  a Tichonov 
condition  (6-58)  for  which  the  kernel  function  defining  the  norm 
/nvfT  is  the  covariance  function,  and  for  which  a = e = 1. 
Thus,  the  relative  weights  of  the  two  summands  (6-58)  are  unique- 
ly determined  in  such  a way  as  to  obtain  an  optimal  (minimum  va- 
riance) result. 

The  observation  equation  to  be  satisfied  is  (6-16): 

1 = AX  + BT  + n . (7-10) 

Putting 


BT  = s'  , (7-11) 

which  is  the  influence  of  the  anomalous  field  on  1,  that  is,  the 
"signal  part"  of  1 , we  have 

1 * AX  + s'  + n , (7-12) 

which  is  the  fundamental  equation  of  least-squares  collocation; 
cf.  eq.(2-l)  of  (Moritz,  1972,  p.7),  with  1 denoted  by  x. 

The  simplest  way  to  derive  (7-8)  and  (7*9)  is  to  use  a 
finite  minimum  norm  principle: 

sTC_1s  + nTD-1n  ■ minimum  , (7-13) 

where  the  vector  s comprises  the  components  of  s'  in  (7-12) 
plus  the  signals  to  be  estimated  In  the  given  problem  (and  only 
these);  C Is  the  covariance  matrix  of  this  vector  s;  cf.  (Mo- 
ritz, 1972,  p . 122) ; a similar  principle  is  used  (ibid. , sec. 2). 


This  approach  Is  completely  equivalent  to  (6-42),  which 
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uses  the  full  Hilbert  space  norm.  Here  only  finite  matrix  opera- 
tions are  required;  Hilbert  space  lurks  invisibly  in  the  back- 
ground (to  use  another  metaphor:  matrix  formulas  such  as  (7-8) 
and  (7-9)  represent  the  finite-dimensional  surface  of  infinitely- 
dimensional  Hilbert  space). 

Thus,  least-squares  collocation  satisfies  two  conditions: 
a minimum  norm  condition  ((6-42)  or  (7-13))  and  a minimum  vari- 
ance condition  ((7-4)  and  (7-5)), 

The  minimum  norm  condition  (especially  in  the  form  (7-13)) 
is  mathematically  convenient  to  handle,  but  physically  it  appears 
somewhat  unnatural:  it  combines  additively  two  different  quanti- 
ties (norms),  which  have  a completely  different  physical  charac- 
ter: the  first  depends  on  the  anomalous  gravitational  field,  the 
second  on  random  measuring  errors.  * 

Much  more  natural  seems  the  minimum  variance  condition: 
the  average  in  (7-3)  (and  similarly  in  (7-5))  is  over  both  signal 
and  noise,  but  affecting  the  same  quantity  (Tp;  - Tp)2:  this  is 
felt  to  be  exactly  as  it  should  be. 

An  interesting  variant  of  the  minimum  variance  approach 
is  to  start  from  a principle  such  as  (7-4)  and  to  postulate  a 
linear  estimation  that  is  symmetric  with  respect  to  the  rotation 
group.  In  this  way  it  is  possible  to  deduce  (7-1),  rather  than 
introducing  it  by  definition  (Sanso.  1978). 

Error  Covariances.  Probably  the  most  important  feature  of 
the  statistical  approach  is  the  possibility  to  estimate,  not  only 
the  values  of  anomalous  field  quantities  and  of  systematic  para- 
meters, but  also  their  accuracy. 

The  error  variances  of  the  parameter  vector  X and  the 
signal  s and  their  covariances  are  expressed  by  the  error  co- 
variance  matrices  Evv,  E and  Ev_,  the  last  being  the  cross- 
covariance  matrix  of  X and  s.  The  vector  s comprises  any 
number  of  signals  (anomalous  field  quantities)  estimated  by  the 
least-squares  collocation  procedure  under  consideration. 

We  have  (Moritz,  1972  , pp. 30-33) 
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xx 

E 

s s 


E 


xs 


(ATr-1A)-1 


c - c E_1c 

ss  sx  xs 

-E  AtHt  , 
xx  ’ 


+ HAE  ATHT 
xx 


where 

H = C TT'1. 

SX 


(7-14) 

(7-15) 

(7-16) 


(7-17) 


The  matrix  Cgg  is  the  signal  covariance  matrix  of  the  s that 
was  to  be  estimated;  E is  its  error  covariance  matrix. 

ss 

These  formulas  express  the  accuracy  of  the  statistically 
optimal  estimates  (7-8)  and  (7-9).  The  statistical  treatment  of 
the  gravity  field,  however,  makes  it  even  possible  to  estimate 
the  accuracy  of  other  linear  estimation  formulas  such  as  (6-61) 
and  (6-62).  We  write  these  equations  in  the  form 


X = G1  , (7-18) 

s = H ( 1 - AX)  = H ( 1 - AG 1 ) 

=■  H ( I - AG ) 1 

- LI  , (7-19) 


I denoting  the  unit  matrix  and  the  vector  s comprising  those 
values  of  T that  were  estimated. 

The  least-squares  estimates  (7-8)  and  (7-9)  may  be  written 
in  an  analogous  form: 


X = G1  , (7-20) 

s = LI  . (7-21) 

Then  equations  (3-52a,b)  of  (Moritz,  1972)  give: 

Exx  * Exx  + - 6WG  - G>T  • (7-22) 

ESS  ‘ ESS  + ‘ * L>T  ‘ (7*23) 


Exx  * Exx  + - 6)Tr<G  - G>T  • (7-22) 

E88  ‘ E88  + ‘ * L)T  ‘ (?*23) 
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Here  E and  E are  the  error  covariance  matrices  of 

XX  ss 

the  estimates  (6-61)  and  (6-62),  and  E and  E are  given 
by  (7-14)  and  (7-15).  It  should  be  noted  that,  in  the  least- 
squares  estimates  (7-8)  and  (7-9),  the  matrices  C and  C 

S X 

are  derived  from  the  covariance  function,  whereas  the  correspond- 
ing matrices  C and  BK  in  (6-61)  and  (6-62)  are  to  be  derived 
from  the  kernel  function  used. 

Geometrical  interpretation  of  minimum  variance . The  con- 
dition of  minimum  variance  has  an  interesting  geometrical  inter- 
pretation by  means  of  the  dual  space,  which  is  the  space  of  all 
bounded  linear  functionals  of  the  elements  of  the  given  linear 
space  (Tscherning,  1978, p. 175). 

Let,  for  simplicity,  be  A = 0 and  D = 0 , that  is,  con- 
sider the  case  of  pure  collocation  without  random  errors  and  sy- 
stematic effects.  Then  (6-6)  reduces  to 

\ = LiT  , ‘ (7-24) 

and  the  least-squares  estimate  (7-9)  becomes 

s = C C'h  . (7-25) 

sx  ' 


Here 


s = FT  (7-26) 

is  a linear  functional,  to  be  estimated,  of  the  anomalous  poten- 
tial T ; as  a matter  of  fact,  (7-25)  could  also  have  been  ob- 
tained by  applying  the  linear  functional  F on  both  sides  of 
eq . (5-13 ) . 

The  estimate  (7-25),  besides  satisfying  the  requirement  of 
minimum  norm  ||T||  , also  satisfies  the  condition 

1 1 F* - F | [ 1 = minimum, 


(7-27) 
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where  ||*IP  denotes  the  norm  in  the  dual  space  and  F*  stands 
for  the  "true"  value  of  the  functional,  of  which  F is  the  esti- 
mate. The  minimum  (7-27)  is  taken  with  respect  to  all  possible 
estimates  of  F as  linear  combinations  of  the  q given  func- 
tionals L . 

In  geometrical  terms,  the  "best"  estimate  (7-25)  is  simply 
the  orthogonal  projection,  of  the  functional  to  be  estimated,  on- 
to the  subspace  of  the  dual  space  spanned  by  the  q given  func- 
tionals L . . 

i 

For  the  case  of  pure  collocation,  this  is  shown  in  (Krarup, 
1978,  p.200).  (The  geometrical  interpretation  is  perhaps  best 
understood  in  the  simple  finite-dimensional  case  of  ordinary  least- 
squares  adjustment,  where  linear  functionals  reduce  to  vectors 
(Moritz , 1966) . ) 

This  geometrical  situation  is  valid  for  arbitrary  kernel 
functions.  A physical  interpretation,  however,  seems  to  be  possible 
only  if  we  identify  the  kernel  function  with  the  covariance  func- 
tion; then  [[F*-  F|[‘  is  simply  the  standard  error  m^,  of  esti- 
mating the  functional  F . In  fact, 

* 

( F* - F)T  = F*T  -FT  = s*-  s = e (7-28) 

F 

is  the  "true  error",  the  difference  between  true  value  s*  and 
predicted  value  s , and 

[|  F*-  F| |* 2 = MU2)  = m2  (7-29) 

is  the  average  of  tp  , that  is,  the  square  of  the  standard  pre- 
diction error  mF  . 

Remark  on  Terminology.  The  minimum  variance  conditions 
(7-4)  and  (7-5)  correspond  to  similar  conditions  in  least-squares 
adjustment  and  least-squares  prediction.  Therefore,  the  name, 
least-squares  collocation,  seems  to  be  appropriate  for  the  case 
in  which  the  solution  is  given  by  (7-8)  and  (7-9). 

In  contrast  to  this  name,  collocation  using  an  arbitrary 
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kernel  function,  which  is  the  case  with  the  methods  described  in 
sections  5 and  6,  might  be  called  collocation  with  kernel  func- 
tions . Pure  collocation  (sec. 5)  would  then  be  kernel  function 
collocation  without  noise  (and  without  systematic  parameters), 
and  the  case  of  "integrated  geodesy" --formul as  such  as  (6-61) 
and  (6-62 ) --woul d then  be  kernel  function  collocation  with  noise, 
including  the  estimation  of  systematic  parameters. 

It  is  evident  that  also  the  results  of  collocation  with  a 
general  kernel  function  satisfy  a minimum  principle  of  form  (6-24) 
and,  dually,  (7-27),  which  might  also  be  called  (geometrical) 
least-squares  principles;  therefore  Krarup  (1978,  p . 1 9 7 ) uses  the 
term  "least-squares  collocation"  also  in  the  case  of  a general 
kernel  function. 

The  terminology  suggested  here  has  the  advantage  of  distin- 
guishing between  "statistical"  and  "geometrical"  collocation. 

We  finally  mention  that  a comprehensive  review  of  quadra- 
tic norm  and  minimum  variance  estimation  principles  is  given  in 
(Grafarend,  1978).  Grafarend  also  considers  rank-deficients  equa- 
tions, occuring  when  the  components  of  the  vector  X are  not  li- 
nearly independent. 
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8 . Review;  Alternatives 


In  this  report  we  have  asked  ourselves  in  which  way  arbi- 
trary geodetic  observations  of  the  earth's  gravitational  field 
can  be  used  to  obtain  information  about  this  field  and  to  detei - 
mine  this  field  to  the  best  possible  extent.  We  have  seen  that 
using  contemporary  mathematical  techniques,  a straight  road  leads 
from  the  nonlinear  observational  equations  to  collocation  with 
kernel  functions  and  least-squares  collocation.  The  main  stages 
on  the  way  were  linearization  (sec. 3)  to  obtain  a linear  impro- 
perly posed  problem  (sec. 4),  to  which  three  different  standard 
methods  of  solution  are  applied:  a restriction  of  the  solution 
space,  leading  to  "pure  collocation"  (sec. 5),  variational  prin- 
ciples of  Tichonov  type,  by  which  measuring  errors  can  be  taken 
into  account,  leading  to  a generalized  collocation  with  kernel 
functions,  and  a statistical  approach,  leading  to  least-squares 
collocation. 

We  shall  now  discuss  these  various  stages  and  possible  al- 
ternati ves . 

L i nearization . All  geodetic  observations  are  nonlinear 
functionals  of  the  potential  V and  of  certain  parameters  X. 
Each  observation  gives  a functional  equation,  and  q observa- 
tions give  a set  of  q of  such  observations.  The  determination 
of  V and  X from  this  system  of  equations  is  a nonlinear  im- 
properly posed  problem.  By  a Taylor  linearization  we  get  a system 
of  q linear  functional  equations  for  the  anomalous  potential 
T and  for  corrections,  again  denoted  by  X,  to  the  parameter 
vector.  For  this  linear  improperly  posed  problem,  standard  mathe- 
matical techniques  exist.  If  the  accuracy  of  linearization  is 
not  sufficient,  we  may  Iterate.  A direct  attack  of  the  nonlinear 
problem  does  not  seem  possible  with  present  mathematical  tools: 
there  is  at  present  no  alternative  to  linearization. 

Why  collocation?  Collocation  In  a mathematical  sense  is 
the  approximation  of  a function  by  fitting  an  analytical  expres- 
sion to  q given  linear  functionals;  for  a very  simple  example 
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see  (Moritz  and  Slinkel  , 1978,  p.32). 

Such  collocation  methods  are  used  in  applied  mathematics 
for  the  approximate  solution  of  differential  equations,  etc. 

(which  are  also  functional  equations!).  Here  the  functionals  are 
usually  supposed  to  be  given  in  a mathematically  exact  way,  and 
the  analytical  expression  is  required  to  fit  these  data  exactly. 

Such  a "pure"  or  "mathematical"  collocation  is,  in  general, 
not  adequately  applicable  to  our  present  geodetic  problem,  in 
view  of  the  inevitable  random  measuring  errors  (noise).  An  excep- 
tion is,  for  instance,  least-squares  interpolation  of  gravity 
anomalies  (interpolation  is  a special  case  of  collocation,  in 
which  the  functionals  are  simply  the  values  of  the  function  at 
discrete  points);  here,  the  measuring  errors  of  gravity  are  con- 
sidered to  be  negligibly  small.  Generally,  measuring  errors,  or 
"noise",  must  be  suitably  taken  into  account:  we  have  a problem 
of  "collocation  with  noise".  (This  is  true  for  the  linearized 
problem,  but  may  be  said  to  hold  even  for  the  original  nonlinear 
problem  because  measurements,  by  their  very  nature,  are  nonlin- 
ear functionals  of  V and  X,  and  the  notion  of  collocation  may 
be  extended,  in  a natural  way,  also  to  the  fitting  of  nonlinear 
functionals.  But,  as  we  have  said,  we  shall  anyway  restrict  our- 
selves here  to  linearized  problems.) 

Thus  the  operational  approach  to  physical  geodesy,  start- 
ing from  the  measurements,  inevitably  lead  to  a collocation  prob- 
lem (with  noise);  there  is,  in  this  sense,  no  alternative  to  col- 
location. What  can  be  chosen  In  different  ways,  is  the  analytical 
expression  for  approximating  the  potential. 

Linear  or  nonlinear  approximation.  Instead  of  a linear 
combination  (5-2)  of  base  functions  one  could,  in  principle,  also 
envisage  other  forms  of  approximation.  For  reasons  of  mathema- 
tical simplicity,  and  also  in  view  of  the  smallness  of  the  quan- 
tities under  consideration  (anomalous  field  quantities  and  cor- 
rections to  the  parameters),  linear  approximations  are  used  prac- 
tically without  exception.  (An  exception  would  be  nonlinear  pre- 
diction to  be  mentioned  later  which  is,  however,  hardly  used  in 
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geodetic  practice.)  From  the  point  of  view  of  theoretical  simpli- 
city and  practical  usefulness,  the  restriction  to  linear  approxi 
mation  methods  seems  to  be  fully  justified. 

Alternative  base  functions.  Besides  kernel  functions,  many 
other  functions  are  being  used  as  base  functions  in  expressions 
such  as  (5-2).  We  only  mention  polynomials  in  polynomial  interpo- 
lation, trigonometric  functions  in  trigonometric  interpolation 
and  approximation  and,  in  geodetic  applications,  multiquadric 
functions  (Hardy,  1976).  Particularly  useful  are  spline  functions 
and  other  finite  elements;  for  an  elementary  discussion  and  com- 
parison with  kernel  functions  cf.  (Moritz,  1978b).  An  interesting 
geodetic  application  of  spline  functions  is  to  a fast  computation 
of  covariance  functions  (Sunkel,  1978). 

These  functions  can  be  very  well  suited  for  particular 
applications,  but  they  cannot  be  used  to  solve  the  general  opera- 
tional problem  of  physical  geodesy  because  they  are  not  harmonic 
outside  the  earth,  which  is  required  for  approximating  the  poten- 
tial T . 

Spherical  harmonic  functions  can  be  used  to  represent  the 
potential,  and  they  are,  in  fact,  fundamental  for  this  purpose. 
From  a practical  point  of  view,  they  are  particularly  suited  for 
representing  the  global  field  at  satellite  altitudes.  The  sample 
functions  of  Giacaglia  and  Lundquist  (1972)  are  finite  linear 
combinations  of  spherical  harmonics  in  a form  convenient  for  cer- 
tain purposes.  For  local  representations  of  the  detailed  gravity 
field,  however,  spherical  harmonics  are  not  applicable.  This  ex- 
cludes their  use  as  base  functions  in  the  present  general  context. 

Kernel  functions  can  equally  well  be  used  for  local  and 
global  purposes;  this  explains  their  application  in  our  general 
operation  approach. 

It  should  be  mentioned  also  that  any  choice  of  base  func- 
tions leads  to  a q x q system  of  linear  equations  for  the  co- 
efficients; for  kernel  functions,  this  system  will  have  a symme- 
tric matrix. 
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Alternative  norms.  Let  us  now  turn  to  the  Tichonov  prin- 
ciple (6-4).  Why  should  we  use  a quadratic  norm,  of  form  (6-19) 
and  (6-21),  respectively?  The  reason  is  simple:  only  then  will 
we  get  a linear  variational  (Euler)  equation,  leading  to  a linear 
combination  of  base  functions. 

This  motivates  the  use  of  a quadratic  norm,  which  is  the 
inner  product  of  T with  itself: 

II  T|f  = (T,T)  , (8-1) 

that  is,  the  use  of  Hilbert  space.  But  why  a Hilbert  space  with 
kernel  functions,  not  some  other  Hilbert  space? 

The  reason  is  that  only  then  the  Euler  equation  will,  in 
fact,  lead  to  a linear  combination  (5-2)  of  base  functions.  Let 
us  illustrate  this  by  means  of  a counterexample.  If  we  restrict 
ourselves  to  the  interpolation  of  functions  f defined  on  a 
sphere  a,  then  a very  obvious  choice  of  norm  would  be 

II  f 1 1 2 = (f.f)  = //f'da  . (8-2) 

a 

The  Hilbert  space  so  defined  does  not  have  a kernel  function  in 
the  proper  tense;  however,  it  may  be  considered  as  a Hilbert 
space  with  a "generalized"  kernel  function,  which  is  a Dirac  del- 
ta function  defined  by 

«(P,Q)  = «(Q,P)  , 

6(P,Q)  =0  if  P 4 Q , (8-3) 

/ / <5  ( P , Q)dop  ■ 1 for  fixed  Q . 

<J 

Then 

(f(P).K(P.Q))  * (f(P),«(P.Q))  - / / f (P)6(P.Q)do  « 

a * 

* f(Q)  / /«(P,Q)dop  « f(Q)  , 
a 
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in  agreement  with  the  definition  (5-19). 

For  our  interpolation  on  the  sphere  we  thus  have,  by  (5-15), 

♦ 3 (P)  - MP.Pj)  • (8-4 ) 

These  "functions"  are  zero  outside  the  interpolation  points.  The 
interpolation  function,  the  linear  combination  (5-2),  has,  there- 
fore, the  property  that  it  is  zero  everywhere  outside  the  points 
at  which  the  functional  values  are  given;  there  is  no  smooth  in- 
terpol ati on . 

This  simple  example  will  illustrate  that  only  quadratic 
norms  with  kernel  functions  can  be  used. 

It  is  not  difficult  to  see  that  a kernel  function  norm 
is  a natural  extension,  to  Hilbert  space,  of  the  usual  quadratic 
norm 

xtK_1x  , (8-5) 

with  a positive  definite  regular  matrix  K,  of  vectors  x in  a 
finite-dimensional  space. 

Thus,  the  Tichonov  approach  with  a quadratic  norm  inevi- 
tably leads  to  collocation  with  kernel  functions;  other  base  func- 
tionscannot  be  obtained  in  this  way. 

In  order  to  further  illustrate  this,  let  us  consider  an- 
other example.  Assume  a function  y * f(x)  to  be  defined  on  the 
unit  circle  (0  = x = Zv),  such  that 

f ( 2ir ) * f(0)  . (8-6) 

« 

Let  us  further  assume  that  the  function  has  a square-integrabl e 
second  derivative  f"(x)  and  define  the  norm  of  f by 

II  ^ II  2 - I2"  |f"(x)|2dx  (8-7) 

o 

(It  Is  actually  a semi-norm,  but  this  Is  Irrelevant  here). 


55 


Expanding  the  function  f(x)  into  a trigonometric  series 

we  have 


1 

f (x)  = aQ  + l (ancosnx  + b^sinnx)  , (8-8) 

' n=l  1 

00 

f"(x)="In2(acosnx+bsinnx)  • (8-9) 

n=  1 n 

The  substitution  of  the  last  equation  into  (8-7)  and  subsequent 
integration  leads  to 

||  f || 2 = tt  l n4(a2  + b2)  . (8-10) 

n=  1 

Using  as  kernel  function 

K(P.Q)  = 7 In_4cos(x  - c)  . (8-11) 

l 

where 


we  have 

,2ir  j2f 

(f(P),  K(P,Q) ) - f 5-1  dx  , (8-12) 

6 dx2  ax2 

which  is  readily  seen  to  be  equal  to  f(c)  , provided 

ao  = 7 /2n  f(x)dx  = 0 • (8-13) 

o 

Thus,  for  functions  satisfying  (8-13),  (8-11)  is  the  reproducing 
kernel.  Interpolation  and  collocation  with  such  a kernel  function 
correspond  to  minimizing  the  norm  (8-7). 

The  minimum  of  a norm  such  as  (8-7),  the  integral  being  ex- 
tended over  an  interval  [a.bj  , is  used  to  define  cubic  spline 
functions  (there  are  nonperiodic  and  periodic  splines);  cf.  (Mo- 
ritz and  Slinkel  , 1978,  pp.26  and  44). 


The  physical  interpretation  of  minimum  norm  (8-7)  is  mini- 
mizing the  total  curvature.  Similarly,  we  could,  in  geodesy,  usn 
the  condition  that  the  total  mean  curvature  of  the  geoid  (the  in- 
tegral, over  the  sphere,  of  the  square  of  mean  curvature)  is  mi- 
nimized, in  order  to  define  a kernel  function. 

Least-squares  collocation  versus  collocation  with  a gene- 
ral kernel  function.  In  least-squares  collocation,  the  kernel 
function  is  chosen  to  be  the  covariance  function;  see  the  end  of 
sec. 7.  the  basic  advantage  of  this  choice  is  the  statistical  sig- 
nificance: it  gives  the  best  linear  estimate  (minimum  variance 
of  the  result). 

Since  the  anomalous  potential  is  not  normally  distributed 
(sec. 7),  the  "best  linear"  estimate  is  not  necessarily  the  abso- 
lutely "best"  estimate:  nonlinear  estimates  may  still  reduce  the 
variance  of  the  result.  Therefore,  nonlinear  prediction  has  been 
considered  (Kaula,  1966;  Grafarend,  1972).  For  most  practical 
purposes,  however,  linear  estimates  seem  to  be  fully  adequate. 

Of  decisive  importance  is  the  possibility,  in  least-squares 
collocation,  to  estimate  the  accuracy  of  the  results.  This  can  be 
done  solely  on  the  basis  of  the  covariances,  without  needing  ac- 
tual measurements.  Therefore,  least-squares  collocation  can  be 
used  for  investigating  possible  configurations  of  measurements, 
for  the  planning  of  surveys,  and  even  for  the  investigation  of 
accuracies  to  be  expected  with  the  use  of  measuring  techniques 
which  are  being  developed  or  only  being  envisaged  for  the  future. 

An  Interesting  feature  has  been  pointed  out  by  Tscherning 
(1977):  if  we  take  the  covariance  function  as  kernel  function, 
then  the  anomalous  potential  T has  an  Infinite  norm,  so  that  T 
Itself  does  not  belong  to  the  Hilbert  space  under  consideration. 

In  physical  terms,  taking  for  example  Interpolation,  this  means 
that  the  interpolating  linear  combinations  of  kernel  functions 
are  all  much  smoother  than  T Itself.  In  practice,  this  Is  an 
advantage  rather  than  a deficiency;  It  Is,  however,  a slight 
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theoretical  drawback  because  it  makes  it  difficult  to  prove  con- 
vergence of  the  approximating  to  the  "true"  function  if  the  den- 
sity of  the  data  approaches  a continuous  distribution  (Tscher- 
ning,  1977).  This  drawback  can  easily  be  remedied  by  changing  the 
higher-order  spherical -harmonic  terms  of  the  covariance  function 
(which  can  anyway  not  be  exactly  determined  because  of  lack  of 
data)  by  an  arbitrarily  small  amount  which  does  not  noticeably 
change  this  function. 

The  exact  covariance  function  could,  of  course,  only  be 
determined  if  we  knew  T exactly  (and  then  we  should  not  need 
collocation!).  So  all  we  can  do  is  to  fit  an  appropriate  analy- 
tical expression  to  empirical  data  (gravity  variances  and  covari- 
ances, degree  variances  from  a spher i cal -harmon i c development, 
etc . ) . 

So,  in  the  practice  of  least-squares  collocation,  we  use 
an  analytical  kernel  function  which  approximates  the  covariance 
function  without  coinciding  with  it.  It  should,  however,  be  close 
enough  so  that  accuracy  estimates  are  meaningful  . 

Even  the  use  of  a general  harmonic  kernel  function  precise- 
ly preserves  the  mathematical  structure  of  the  gravity  field,  ex- 
pressed by  relations  between  the  potential,  gravity  anomalies,  de- 
flections of  the  vertical,  etc.. 

The  approach  followed  in  practice  by  most  workers  in  this 
field,  whether  they  explicitly  favor  a statistical  i nterpretation 
or  not,  is  to  use  a kernel  function  which  has  a simple  analytical 
expression  and  exhibits  the  main  features  of  an  empirical  covari- 
ance function. 

Relation  to  least-squares  adjustment.  Least-squares  collo- 
cation differs  from  least-squares  adjustment  in  two  respects. 

First,  from  a physical  point  of  view,  adjustment  contains 
only  one  kind  of  statistical  quantity,  namely  the  random  measur- 
ing errors  (noise);  in  least-squares  collocation  there  are  two 
physically  different  quantities  that  are  treated  statistically; 
quantities  of  the  anomalous  gravity  field  (the  signal)  and  measur- 
ing noise. 
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Second,  from  a formal -mathemati cal  point  of  view,  the  im- 
portant distinction  is  that  the  signal,  in  contrast  to  the  noise, 

is  a continuous  function;  therefore,  an  adequate  treatment  requ i r- 

> 

es  some  i nf i ni tely-dimensional  function  space,  especially  Hilbert 
space.  Operations  in  Hilbert  space  show  formal  similarities  with 
operations  in  finite-dimensional  space  but  are  qualitatively  dif- 
ferent, just  as  differential  equations  are  qualitatively  different 
from  difference  equations  although  there  are  formal  similarities. 

It  is  true  that  the  gravity  field  at  satellite  elevations 
can  be  adequately  described  by  a spherical -harmonic  expansion  tru.  - 
cated  at  a sufficiently  high  degree;  the  coefficients  of  such  a de- 
velopment do  form  a finite-dimensional  vector.  Thus,  if  we  exclu- 
sively work  at  satellite  altitudes,  we  might,  in  fact,  replace  Hil- 
bert space  by  a finite-dimensional  space,  without  essentially  im- 
pairing the  accuracy. 

This  situation  changes  essentially  if  we  consider  the  gra- 
vity field  at  the  earth's  surface  by  including  terrestrial  obser- 
vations or,  e.g.,  satellite  altimetry.  The  detailed  gravity  field 
at  the  earth's  surface  cannot  be  adequately  described  by  a spher- 
ical-harmonic expansion,  neither  from  a theoretical  point  of  view 
--because  the  convergence  cannot  be  guaranteed--nor  from  a practi- 
cal point  of  view--because,  if  at  all  possible,  such  an  expansion 
would  require  an  excessively  high  number  of  terms,  which  is  be- 
yond the  capacity  of  any  present  digital  computer. 

Hence,  the  general  replacement  of  Hilbert  space  by  a finite- 
dimensional space  is  neither  theoretically  nor  practically  feasible. 

It  is  also  not  necessary  since,  as  we  have  seen,  the  practi- 
cal collocation  formulas  are  finite-dimensional  matrix  formulas. 

The  approach  of  (Moritz,  1972)  works  entirely  with  finite  matrices; 
the  Hilbert  space  character  expresses  Itself  only  In  the  fact  that 
covariances  are  propagated,  not  by  matrix  operations,  but  by  linear 
operations  (such  as  differentiation)  of  covariance  functions.  See 
also  (Moritz,  1976,  sec. 4). 

A feature  of  collocation  that  Is  theoretically  most  Inter- 
esting is  that,  from  a formal -geometri cal  point  of  view,  least- 
squares  collocation  can  be  considered  as  a problem  of  least-squares 
adjustment  In  Hilbert  space  (Krarup,  1969,  pp. 34-41). 
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We  also  mention  that  the  usual  least-squares  solution  of 
overdetermined  systems  of  linear  equations  can  be  considered  as 
a "quasi-solution"  according  to  the  theory  of  improperly  posed 
problems  (Tikhonov  and  Arsenin,  1977,  Chapter  I). 

Discrete  versus  continuous  data.  In  general,  geodetic  data 
are  discrete  measurements  of  a finite  number;  the  present  report 
is  based  on  this  situation.  There  are  exceptions,  for  instance, 
continuous  recordings  of  data  profiles;  they  are  usually  convert- 
ed into  discrete  data  by  a representative  selection. 

A practically  more  important  case  is,  for  instance,  a very 
dense  gravity  survey  around  a station  at  which  precise  deflections 
of  the  vertical  are  to  be  computed.  In  such  cases,  a combination 
of  collocation  with  other  methods,  such  as  integral  formulas  (Mo- 
ritz, 1975),  might  be  appropriate. 
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